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Abstract 

A  numerical  investigation  of  the  transport  phenomena  and  performance  of  a  proton  exchange  membrane  fuel  cell  (PEMFC)  with  various  design 
parameters  of  the  transport  component  is  presented.  A  three-dimensional  fuel  cell  model,  incorporating  conservations  of  species,  momentum,  as 
well  as  current  transport,  is  used.  The  Bulter-Volmer  equation  that  describes  the  electrochemical  reaction  in  the  catalyst  layer  is  introduced;  the 
activation  overpotential  connects  the  solid  phase  potential  field  to  that  of  the  electrolyte  phase.  Through  cell  performance  simulation  with  various 
channel  aspect  ratios  and  gas  diffusion  layer  (GDL)  thicknesses,  a  slender  channel  is  found  suitable  for  cells  operating  at  moderate  reaction  rate, 
and  a  flat  channel  produces  more  current  at  low  cell  voltage.  Plots  of  transverse  oxygen  concentration  and  phase  potential  variation  indicate  that 
these  oppositely  affect  the  local  current  density  pattern.  The  relative  strengths  of  these  two  factors  depend  on  the  transport  component  position  and 
geometry,  as  well  as  on  the  cell  operating  conditions.  Consequently,  the  curves  of  cell  output  current  density  demonstrate  that  the  optimal  GDL 
thickness  increases  as  the  cell  voltage  decreases.  However,  at  the  lowest  considered  cell  voltage  of  0.14  V,  optimal  thickness  decreases  as  that  of  a 
thick  GDL.  The  oxygen  deficiency  caused  by  long  traveling  length  and  clogging  effect  of  liquid  water  reverses  this  relationship. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Fuel  cell  technology  greatly  impacts  industry  development 
and  has  attracted  recent  attention  because  the  shortage  of  fos¬ 
sil  oil  and  environmental  protection  have  become  increasingly 
important  issues.  The  proton  exchange  membrane  fuel  cell 
(PEMFC)  particularly  has  excellent  high  efficiency  and  zero 
emission  characteristics  and  can  be  operated  at  low  temperature. 
It  is  a  commonly  proposed  substitute,  therefore,  to  the  internal 
combustion  engine  [1]. 

The  PEMFC  uses  pressurized  reactants  that  circulate  through 
grooved  bipolar  plate  channels  to  provide  cell  energy.  Oxidation 
occurs  on  the  anode  side,  releasing  energetic  electrons  and  pro¬ 
tons.  Meanwhile,  the  cathode  reduces  oxygen  with  electrons  and 
protons  to  form  water.  Electric  energy  is  released  during  elec¬ 
trons  transport  through  external  loading.  Structural  simplicity, 
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compact  operating  principles,  and  product  neatness  are  the  most 
attractive  features. 

PEMFC  performance  is  dictated  by  several  factors  such  as 
transport  component  geometry  and  morphology,  and  operating 
conditions.  Giorgi  et  al.  [2]  investigated  diffusion  layer  porosity 
influence  on  effective  catalyst  activity  for  the  cathode  reduction 
reaction.  Thickness  and  diffusion  layer  structure  effects  on  low 
Pt-loading  electrode  performance  for  PEMFC  were  also  con¬ 
ducted  experimentally.  An  optimal  thickness  of  diffusion  layer 
was  found  due  to  its  lower  electrical  resistance.  Decreased  per¬ 
formance  with  the  thickest  diffusion  layer  was  attributed  to  a 
long  reactant  transport  passage  and  to  the  flooding  problem. 
Jordan  et  al.  [3]  reported  the  carbon  support  and  electrode  load¬ 
ing  effects  on  PEMFC  performance.  Various  parameter  impacts, 
such  as  diffusion  layer  thickness,  PTFE  content  and  morphology, 
as  well  as  operating  temperature  and  humidification  condition 
on  cell  performance  were  addressed.  Lee  et  al.  [4]  asserted  from 
experimental  results  that  an  optimal  bolt  torque  was  obtained 
for  a  soft  commercial  diffusion  layer  because  of  porosity  and 
electrical  contact  resistance  changes.  Their  results  also  revealed 
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molar  concentration  (molm-3) 
form  drag  constant 
diffusivity  (m2  s-1) 

Faraday  constant,  96,500  C  mol-1 
gravitational  acceleration  (ms-2) 
height  (m) 

current  density  (A  m-2) 
relative  permeability 
permeability  (m2) 
length  (m) 

catalyst  loading  (mg  m-2) 
molecular  weight  (kg  mol-1) 
number  of  species 
pressure  (Pa) 

number  of  electron  transferred 
rate  constant 

universal  constant,  8.314  J  mol-1  K-1 

saturation  level 

source  term 

temperature  (K) 

mixture  velocity  (ms-1) 

pore  volume  of  porous  medium  (m3) 

Voltage  (V) 

mass  fraction  of  species 
width  of  domain  (m) 
molar  fraction 
switching  function 


Greek  symbols 
p  transfer  coefficient 

y  water  content  in  membrane 

8  tortuosity 

£  dry  porosity 

r\  overpotential  (V) 

0  contact  angle 

l±  dynamic  viscosity  (N  s  m-2) 

v  kinematic  viscosity  (m2  s-1) 

§  stoichiometry  of  hydrogen  or  oxygen 

a  electric  conductivity  (S  m-1) 

q  surface  tension  (Nm-1) 

0  potential  (V) 

cp  concentration  dependence 


Superscripts 
c  capillary 

eff  effective  value 


Subscripts 
a  anodic 

c  cathodic 

cat  catalyst 

cel  cell 

cha  channel 

con  condensation 


eva 

evaporation 

F 

Forchheimer  term 

g 

gaseous  phase 

j 

volumetric  transfer  current  density 

k 

anode  or  cathode 

1 

liquid  phase 

m 

momentum 

mas 

per  unit  mass  of  carbon  supported  catalyst  particle 

mem 

membrane  phase 

P 

index  for  phases 

pit 

platinum 

s 

species 

sat 

saturation 

sol 

solid  phase 

V 

per  unit  volume 

w 

water 

a 

index  for  species 

that  a  cell  with  high  bore  torque  corresponding  to  a  thin  diffusion 
layer  generates  more  current  at  a  moderate  reaction  rate,  and  a 
moderate  bore  torque  performs  better  at  a  high  reaction  rate  for 
this  kind  of  diffusion  layer. 

Several  elaborate  mathematical  models  were  developed  to 
provide  qualitative  insights  into  the  transport  phenomena  in  a 
PEMFC,  besides  the  experimental  works.  He  et  al.  [5]  pre¬ 
sented  a  two-dimensional,  two-phase,  multi-component  trans¬ 
port  model  investigating  gas  and  liquid  water  hydrodynamic 
effects  on  PEMFC  performance  employing  an  inter-digitated 
flow  field.  Parametric  study  on  electrode  thickness  effect  indi¬ 
cated  that  within  an  electrode  finite  thickness,  average  cur¬ 
rent  density  increases  monotonically  with  electrode  thickness 
increase,  no  matter  what  the  operating  voltage.  Furthermore, 
beyond  this  finite  thickness,  cell  performance  decreases.  This 
phenomenon  was  explained  by  reactant  transport  and  liquid 
water  removal  viewpoints.  The  effects  of  channel  and  shoulder 
width  ratio  (CIS  ratio)  revealed  that  cell  performance  enhances 
with  a  greater  CIS  ratio.  Natarajan  and  Nguyen  [6]  proposed 
a  half-cell,  transient  model  for  the  cathode  using  the  conven¬ 
tional  flow  field  to  evaluate  design  parameter  influences  on  fuel 
cell  transport  processes.  Contrary  to  He  et  al.  [5],  they  found, 
as  diffusion  layer  thickness  decreases  at  lower  cathode  overpo¬ 
tential,  the  reactant  transport  is  not  affected  by  the  liquid  water 
and  at  some  diffusion  layer  thickness,  current  density  loss  at  the 
shoulder  region  outweighs  channel  region  gain,  leading  to  an 
optimum  diffusion  layer  thickness.  However,  a  thinner  diffusion 
layer  thickness  results  in  better  performance  and  an  optimum 
thickness  does  not  exist  at  a  higher  reaction  rate. 

In  the  two-dimensional  half-cell  model  performed  by  West 
and  Fuller  [7],  effects  of  rib  sizing  and  geometry  of  GDL  on  cur¬ 
rent  and  water  distributions  within  the  cell  were  reported.  Special 
attention  was  given  to  the  relation  between  rib  width  and  water 
content  in  the  membrane.  Kulikovsky  et  al.  [8]  conducted  the 
parametric  studies  with  high  and  low  values  of  solid  phase  elec¬ 
trode  conductivity  elucidating  the  possibility  of  locally  reducing 
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catalyst  loading  in  the  catalyst  layer.  Their  study  concluded  that 
at  low  solid  phase  conductivity,  cell  reaction  concentrates  mainly 
at  the  rib  region  catalyst  layer  and  the  catalyst  can  be  removed 
from  the  channel  region  catalyst  layer  with  little  impact  on  total 
cell  performance.  Simplified  descriptions  of  the  kinetic  and  mass 
transport  model  developed  from  Jeng  et  al.  [9]  concluded  that 
GDL  effectiveness  decreases  with  reaction  rate  and  increases 
with  flow  channel  width.  Also,  when  GDL  porosity  is  low,  cell 
performance  decreases  with  an  increase  in  GDL  thickness.  Sun 
et  al.  [10,11]  conducted  a  series  of  investigations  on  catalyst 
layer  structural  parameter  influences  with  an  improved  cross- 
the-channel  model.  Various  design  and  operating  factors  such  as 
C/S  ratio,  orthotropic  GDL  conductivity,  and  electrode  compres¬ 
sion  effects  were  investigated  in  the  model  while  charge  transfer 
and  oxygen  diffusion  were  simultaneously  taken  into  account. 
Gas  and  electron  transport  plays  an  essential  role  in  cathode 
performance  in  that  both  electric  conductivity  and  GDL  thick¬ 
ness  could  be  vital  parameters  for  transport  component  optimal 
design,  as  conclusions  indicated. 

Most  groups  researching  cell  performance  simulation  have 
recently  begun  applying  commercial  or  self-developed  compu¬ 
tational  fluid  dynamic  (CFD)  algorithms  in  their  study.  Hence,  a 
more  rigorous  three-dimensional  simulation  has  been  made  pos¬ 
sible  for  investigating  the  complex  transport  phenomena.  The 
species  and  properties  in  the  physical  domain  can  be  described 
in  greater  detail,  and  numerical  simulation  results  are  more  reli¬ 
able  than  one-  or  two-dimensional  models  [12-17]. 

The  objective  of  this  study  is  to  present  the  interactions 
between  transport  component  design  and  PEMFC  performance. 
A  three-dimensional,  multi-species  and  two-phase  model  is  used 
to  investigate  channel  aspect  ratio  and  GDL  thickness  effects  on 
cell  performance.  Local  resolution  of  current  generation  domi¬ 
nant  mechanisms  is  also  emphasized,  revealing  electrochemical 
crucial  reaction  factors,  such  as  reactant  concentration  and  acti¬ 
vation  overpotential. 

2.  Mathematical  formulations 

Fuel  cell  configuration  is  symmetric  with  a  typical  straight 
flow  field,  so  only  half  of  the  domain  is  considered  here.  The 
model  considers  nine  components,  as  Fig.  1  illustrates.  The 
anode  and  cathode  catalyst  layers  (CLs)  and  membrane  are  con¬ 
sidered  separate  entities,  despite  the  ultra- thinness  of  the  CLs. 
Two  GDLs  are  compressed  by  current  collector  ribs  on  each  cell 
side  to  deliver  reactants  to  the  reaction  sites  under  the  ribs.  The 
spaces  in  the  bipolar  plate  (BP)  grooves,  commonly  called  the 
anode  and  cathode  channels,  are  passages  in  which  gases  are 
transported  throughout  the  cell.  This  work  describes  the  follow¬ 
ing  transport  phenomena: 

•  three-dimensional  convection  and  diffusion  in  flow  channels 

and  porous  media; 

•  electrochemical  reactions  in  the  catalyst  layers  of  the  anode 

and  the  cathode; 

•  multi-species  mass  transport; 

•  formation  and  transport  of  liquid  water; 


Cathode  BP 

Cathode 

Channel 

Cathode  GDL 

Cathode  CL 
Membrane 

Anode  CL 
Anode  GDL 


Anode 

Channel 

Z 


Fig.  1.  Physical  domain  of  interest,  the  channel  height  is  in  Z-direction  and  the 
channel  width  is  in  F-direction. 


•  electronic  current  transport  through  bipolar  plate  and  elec¬ 
trode; 

•  positive  ionic  current  transport  through  the  membrane  and 
catalyst  layer. 

2.7.  Basic  assumption 

Basic  assumptions  are  made  to  simplify  actual  cell  condi¬ 
tions  in  the  theoretical  model  and  thus  facilitate  the  modeling 
approach  of  transport  component  influence  on  transport  phe¬ 
nomena  and  cell  performance.  The  following  are  the  most  impor¬ 
tant: 

•  the  Reynolds  number  of  the  fluid  is  below  100  due  to  low 
mixture  velocity,  laminar  flow  is  considered; 

•  the  GDLs,  CLs  and  membrane  are  modeled  as  porous  media; 

•  each  component  has  isotropic  transport  properties; 

•  the  gaseous  phase  of  the  working  fluid  behaves  as  an  ideal 
gas; 

•  the  electric  potential  on  the  outer  surface  of  each  bipolar  plate 
is  constant; 

•  no  charge  accumulates  in  the  electrodes  and  the  domain  is 
electrically  neutral; 

•  the  system  operates  in  a  steady  and  isothermal  state. 

2.2.  Governing  equations 

A  unified  domain  approach  is  employed  to  avoid  tedious 
boundary  condition  appointments  on  the  component  interface  in 
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the  following  formulations.  Instead,  proper  transport  properties 
such  as  porosity  and  permeability  are  designated  on  each  dis¬ 
tinct  domain.  This  method  can  significantly  simplify  the  model 
system  equations  and  numerical  processes. 

The  continuity  equation  is  used  to  describe  the  mass  conser¬ 
vation  of  mixture  throughout  the  domain: 


V(sefipU)  =  0 


(1) 


where  the  mixture  density  p  is  the  volume- weighted  average  of 
the  phase  mass  concentration  for  the  consideration  of  the  two- 
phase  flow  [18].  It  can  be  expressed  as 


p  =  p\s  T-  pg(l  -  s ) 


(2) 


where  s  stands  for  the  saturation  level  in  porous  medium,  repre¬ 
senting  the  volume  fraction  of  the  pore  occupied  by  the  liquid 
phase  and  can  be  given  as 


.v  = 


v\ 


V\  +  vg 


(3) 


The  generalized  Navier-Stokes  equation  is  introduced  to  rep¬ 
resent  momentum  conservation  of  the  mixture;  a  source  term  is 
included  to  consider  the  additional  drag  forces  in  the  porous 
medium: 


,eff , 


eff. 


,eff  ..effi 


,eff 


£eff  =  £(1  —  s),  jueff  is  the  effective  viscosity  of  the  mixture,  pk 
the  kinetic  density,  Sm  is  the  sum  of  Darcy  and  Forchheimer  drag 
forces  associated  with  the  morphology  (porosity  and  permeabil¬ 
ity)  of  the  porous  medium: 


eff 


Pls  +  Pg(l  —  s) 

(5) 

{kr\/V\)  +  (fcrg/Vg) 

V  V 

(6) 

P\kr\  T  Pgkrg 

VI  Vg 

(Sefi)2U  ( Meff  eff  \ 

-  4k  U  +  £  Wl 

(7) 

The  momentum  equation  automatically  becomes  the  Darcy 
equation  in  the  porous  area  where  the  magnitudes  of  convection 
and  diffusion  term  are  extremely  small.  In  Eq.  (5),  the  relative 
permeabilities  of  the  liquid  and  gas  phase  can  be  expressed  as  a 
function  of  saturation  [18]: 


kr  i  =  s3 
kre  =  (1  -  S)3 


(8) 

(9) 


The  mass  fraction  of  each  species  in  the  gas  mixture  can  be 
given  as  follows: 


V  I  £eff pu>aUkrg  —  )  =  V [p A* (A11  )"Vwa]  +  S: 


V 


eff  A- 


(10) 


v 


Notably,  the  diffusivity  is  modified  from  that  of  the  free 
stream  value  by  applying  the  Bruggemann  correction.  The 
parameter  8  is  the  tortuosity  of  the  porous  medium.  The  source 


term  Ss  defines  the  creation  or  consumption  of  species  at  the 
electrode  catalyst  sites  and  is  given  by 


So,  = 


Sj,  cMq: 


Sr2o  = 


4000 F 

—Sj,cMu2o 
2000 F 


+  Si 


Svh  = 


Sj,&Mur 


2000 F 


(11) 


(12) 


(13) 


At  the  electrode  catalyst  sites,  reactants  undergo  an  electro¬ 
chemical  reaction.  Hydrogen  is  oxidized  and  oxygen  reduced 
at  the  anode  and  the  cathode,  respectively.  These  two  reactions 
are  driven  by  the  potential  difference  between  the  solid  phase 
and  the  electrolyte  phase,  called  the  activation  overpotential  77. 
The  activation  overpotential  in  the  anode  tends  to  energize  the 
electrons  lost  by  hydrogen;  at  the  cathode,  hydrogen  ions  and 
electrons  react  with  oxygen  to  form  water.  The  Bulter-Volmer 
equation  describes  this  important  phenomenon  and  is  related  to 
the  source  terms  in  Eqs.  (1 1)— (13): 


Sj,k  —  Av?"  ref.  A; 


exp 


c 


( Pa 


a 


(4) 

X 

by 

where  Av 

Qy,ref 

Pk,SiFr]k 

RT 


exp 


-Pk,cFr]k 

RT 


(14) 


particle  and  is  closely  related  to  the  morphology  and  platinum 
loading  of  the  catalyst  layer.  An  expression  for  the  relationship 
among  these  parameters  is  written  as 


Av  = 


A  mas  ^  pit 


Ht 


(15) 


cat 


where  Amas  is  the  catalyst  area  per  unit  mass  of  catalyst  parti¬ 
cle,  mpit  represents  the  catalyst  loading  and  Hcat  stands  for  the 
catalyst  layer  thickness. 

The  term  /ref  ^  in  Eq.  (14)  is  the  exchange  current  density, 
which  characterizes  the  catalyst  layers.  Due  to  sluggishness 
of  oxygen  reduction  reaction  at  the  cathode,  its  value  is  sev¬ 
eral  orders  of  magnitude  smaller  than  that  at  the  anode.  Hence, 
the  cathode  exhibits  significant  activation  overpotential  during 
cell  operation.  The  experimental  results  of  Parthasarathy  et  al. 
yielded  the  following  relation  [19]: 


.  ,  4001  \ 

*ref,02  =  10  exP  (  3.507  - —  ) 


(16) 


Two  charged  species,  electrons  and  protons,  are  transported 
in  the  fuel  cell  and  subjected  to  individual  driving  forces  deter¬ 
mined  by  the  potential  gradient.  The  negatively  charged  elec¬ 
trons  flow  through  the  catalyst  layer,  the  gas  diffusion  layer, 
and  the  bipolar  plate;  the  positively  charged  protons  flow  goes 
through  the  catalyst  layer  and  the  membrane.  The  catalyst  layer 
is  a  source  of  charges  at  the  anode  and  a  sink  at  the  cathode.  The 
assumption  of  electro-neutrality  yields  the  following  expres¬ 
sions  for  electrical  current  conservation  in  the  catalyst  layer: 


^  (  <TsolV0sol)  —  Sj,k 


(17) 
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^mcm^7  0rncm)  —  $j,k  (18) 

Outside  the  catalyst  layer,  no  current  sink  or  source  is  pre¬ 
sented  and  the  right-hand  sides  of  Eqs.  (17)  and  (18)  are  equal 
to  zero,  suggesting  that  no  species  is  created  or  consumed;  elec¬ 
trical  current  conservation  is  therefore  easily  understood. 

In  the  membrane,  the  ionic  conductivity  <7mem  is  strongly 
related  to  the  water  content  y,  defined  as  the  ratio  of  the  number 
of  water  molecules  to  the  number  of  the  charge  sites  [20] : 

(Tmem  =  (0.005139y  —  0.00326)  | exp 

(19) 

An  empirical  relationship  between  the  water  content  in  the 
membrane  and  the  partial  pressure  of  the  water  is 


coefficient.  The  first  term  on  the  right-hand  side  describes  the 
capillary  pressure  effect,  whereas  the  second  term  is  the  gravi¬ 
tational  separation  between  the  liquid  and  gas  phases  of  water. 
The  variable  Dc  is  the  capillary  diffusion  coefficient  and  can  be 
expressed  as  a  function  of  saturation  level,  surface  tension,  and 
permeability  [22,23]: 


/  peff  \  1/2  9 

JS  kp\kf  g 

cos 0c(£k)  (  3.789s2+3.338s  0.966) 

zr=- 


kr  iVg  +  kr  gVi 


(24) 


where  kr i,  krg  are  the  relative  permeabilities  associated  with  pore 
space  reduction  by  the  co-existence  of  multiphase  fluids. 

The  last  term  in  Eq.  (23)  is  introduced  to  account  for  liquid 
water  condensation  or  evaporation,  and  can  be  expressed  as  [5]: 


0.043  +  17.81a  —  39.85a2  +  36a3,  for  0  <  a  <  1 
12.6  +  1.4a,  for  1  <  a  <  3 

(20) 


where  a  is  the  water  activity  and  is  given  by 

*wP 
a  =  - 

Psat 


(21) 


In  the  above  equation,  the  saturation  pressure  varies  with 
temperature  and  can  be  determined  from  the  thermodynamic 
table  or  using  the  following  empirical  expression: 

logic  Psat  =  -2.1794  +  0.02953T-  9.1837  x  10“5r2 

+  1.4454  x  10“7r3  (22) 


where  T  is  in  unit  of  °C  and  psat  is  in  unit  of  bars.  This  formulation 
can  be  utilized  with  a  relative  error  less  than  1.5%  in  the  range 
of  20-100  °C. 

During  fuel  cell  operation,  water  partial  pressure  in  the  elec¬ 
trode  may  exceed  its  saturation  pressure  if  the  local  water  con¬ 
centration  is  high.  Therefore,  liquid  water  is  possibly  formed 
and  occupies  the  electrode  pore  space.  Operating  the  cell  at 
a  high  reaction  rate  may  cause  severe  mass  transport  overpo¬ 
tential  because  the  diffusion  species  are  blocked.  Additionally, 
extremely  small  pores  in  the  fuel  cell  porous  media  cause  the  cap¬ 
illary  force  to  dominate  liquid  water  transport  in  the  hydrophilic 
surfaces.  Smaller  pores  correspond  to  stronger  capillary  forces. 
However,  the  actual  expression  of  this  force  cannot  be  formu¬ 
lated  because  the  real  liquid-gas  interface  configuration  is  not 
tractable.  The  generalized  Richards  equation,  originally  devel¬ 
oped  for  the  oil  recovery  field  and  adapted  by  Wang  and  Becker- 
mann  [21]  to  study  two-phase  flow  transport  in  capillary  porous 
medium,  is  applied: 

V  (stnpUkrX/  \  +  VM 

=  V  (eeffDcVs  -  P/lSlSl  ~  pg)g)  +  $,  (23) 

\  kp\vg  A  kpgv\  J 


£effXw 

S\  —  4#[  rcon  (xwp  —  PsciOy 

K1 

rr 

+  reva£e  Sp\ (xw P  ~  PsatXl  -  y) 


(25) 


where  y  is  a  constant  with  unity  or  zero  value  depending  on  water 
species  condensation  or  evaporation  scenarios  in  the  porous 
medium. 


2.3.  Boundary  conditions 


Model  boundaries  fall  into  three  categories  as  indicated  by 
Fig.  1 — symmetric  boundaries  (SBs),  impermeable  boundaries 
(IBs),  and  channel  boundaries  (CBs).  The  SB  conditions  are 
imposed  on  both  transverse  sides  of  the  model  and  physical 
quantities  such  as  mass  flux  or  momentum  flux  have  zero  gra¬ 
dients.  The  Neumann  conditions  are  also  assigned  on  IBs,  the 
top  (or  bottom)  and  side  surfaces  of  the  channels  adjacent  to  the 
bipolar  plates.  Additionally,  non-slip  conditions  are  applied  for 
the  IB  velocity  fields. 

Unlike  these  simple  boundary  conditions,  the  CB  conditions 
importantly  determine  operating  conditions  and  cell  perfor¬ 
mance.  At  the  inlet  ports  of  the  channels,  velocities,  thermo¬ 
dynamic  states,  and  mass  fractions  of  mixture  are  specified 
according  to  the  desired  operating  scenarios  and  mixture  of  inter¬ 
est.  Based  on  the  stoichiometric  flow  ratio,  the  inlet  velocity  of 
the  reactant  is 


§  IW cq\L  cha  RTgm 
q FW cha^cha  (l  -  XaMPk 


(26) 


The  mass  fractions  of  the  species  on  the  CBs  are  specified 
according  to  operating  pressure  and  fully  saturated  humidifi¬ 
cation  condition.  Dalton’s  law  and  the  ideal  gas  law  yield  the 
molar  fraction  of  individual  species.  The  mass  fraction,  used  in 
the  preceding  formulation,  is  given  by 


(27) 


where  N\  represents  the  liquid  water  flux  due  to  electro-osmosis 
in  the  membrane,  it  is  associated  with  proton  transport  and  is 
a  function  of  local  current  density  and  the  electro-osmotic  drag 


Specific  solid  phase  potential  is  assigned  to  each  bipolar 
plate  outer  surface  in  addition  to  these  transport  variables.  This 
value  is  generally  set  to  zero  on  the  anode  BP  outer  surface, 
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Table  1 


Cell  geometries  and  operating  conditions  of  base  model 


Geometry  and  condition 

Value 

Unit 

Domain  length 

50 

mm 

Domain  width  (half) 

0.8 

mm 

Gas  channel  width  (half) 

0.4 

mm 

Domain  height 

4.303 

mm 

Gas  channel  height 

0.8 

mm 

Diffusion  layer  height 

0.254 

mm 

Catalyst  layer  height 

0.01 

mm 

Membrane  height 

0.175 

mm 

Molar  ratio  of  cathode  side  dry  air  (N2/O2) 

79/21 

Reactant  relative  humidity 

100% 

Reactant  stoichiometry 

3 

Cell  back  pressure 

2 

atm 

Operating  temperature 

353 

K 

and  a  cell  total  overpotential  is  set  on  the  cathode  BP.  Mean¬ 
while,  a  membrane  phase  potential  zero  gradient  is  applied  at  the 
interface  between  GDL  and  CL,  representing  a  protonic  current 
absence  through  this  interface  because  of  the  lack  of  a  conduct¬ 
ing  medium.  This  method  yields  the  local  current  density  in  the 
cell  from  ohmic  law;  it  also  accurately  yields  the  catalyst  layer 
activation  overpotential.  The  cell  potential  can  be  obtained  from 
the  following  expression: 

Teel  —  Voc  T  Otot  (28) 

where  rjto t  is  the  total  cell  overpotential  and  Voc  is  the  open  circuit 
voltage  given  by  [17]: 

Voc  =  0.2329  +  0.00257  (29) 

2.4.  Method  of  solution 

The  nonlinearity  and  coupling  in  the  model  formulation 
eliminate  an  analytical  solution  possibility.  The  forgoing  trans¬ 
port  equations  are  converted  into  a  conservation  form  of 
convection-diffusion  conservation  equations  and  numerically 
solved  using  the  commercial  computational  fluid  dynamic  code 
CFD-ACE+,  based  on  the  control  volume  formulation  and  the 
SIMPLE  algorithm  [24-28]. 

Three  structure  mesh  systems — 60  x  30  x  12,  80  x  40  x  16, 
and  90  x  48  x  18  are  constructed  to  explore  numerical  result 
dependence  on  computational  cell  numbers.  Base  model  geome¬ 
tries  and  operating  conditions  are  listed  in  Table  1.  Model 


Table  2 


Electrochemical  parameters  and  transport  properties 


Parameters  and  properties 

Value 

Unit 

Sources 

Porosity  of  the  diffusion  and  catalyst 

0.4 

[29] 

layer 

Porosity  of  the  membrane 

0.28 

[29] 

Permeability  of  the  diffusion  and 

2.3E-11 

2 

m 

[26] 

catalyst  layer 

Permeability  of  the  membrane 

1.0E-18 

2 

m 

[26] 

Tortuosity  of  the  diffusion  and 

1.5 

[30] 

catalyst  layer 

Tortuosity  of  the  membrane 

3 

[30] 

Condensation  rate  constant 

100 

s-1 

[5] 

Evaporation  rate  constant 

100 

atm-1  s-1 

[5] 

Concentration  dependence  of  H2 

0.5 

[31] 

Concentration  dependence  of  O2 

1 

[31] 

Transfer  coefficients  at  anode  (anodic 

0.5 

[30,32] 

and  cathodic) 

Transfer  coefficients  at  cathode 

1.5 

[30,32] 

(anodic  and  cathodic) 

Electrical  conductivity  of  electrode 

114 

Sm  1 

[33] 

Catalyst  loading 

0.4 

mg  cm-2 

Catalyst  surface  area  per  unit  mass 

100 

2  -l 

m  g 

[34] 

Contact  angel 

0 

O 

[18] 

Exchange  current  density  for  anode 

1.4E5 

Acm~3 

[35] 

reaction 

Reference  concentration  of  oxygen 

3.39E-6 

mol  cm  3 

[35] 

Reference  concentration  of  hydrogen 

5.64E-5 

mol  cm-3 

[35] 

component  electrochemical  parameters  and  transport  proper¬ 
ties,  obtained  mainly  from  open  literature  [29-35],  are  listed  in 
Table  2.  Numerical  calculations  are  performed  on  a  Pentium  IV 
2.4  GHz  PC  with  a  1G  RAM.  The  mesh  system  (80  x  40  x  16)  is 
adopted  because  current  density  values  are  satisfactory  with  an 
error  range  of  2%  when  using  base  model  geometries  and  param¬ 
eters  of  these  two  tables,  see  Table  3.  Based  on  the  reported 
operating  conditions  and  component  geometries  of  Um  and 
Wang  [17],  data  in  Table  4  show  that  the  present  model  to  be  con¬ 
sistent  with  previous  studies.  However,  as  pointed  by  Sivertsen 
and  Djilali  [36],  this  is  not  sufficient  to  validate  computational 
model,  comparison  with  experimental  data  providing  detailed 
information  of  key  quantities  is  essential  in  the  future. 

3.  Results  and  discussion 

Transport  component  design  effects  on  physical  quantity  dis¬ 
tributions  and  PEMFC  performance  are  analyzed  by  numerical 


Table  3 

Results  of  three  computation  grid  systems  based  on  the  parameters  in  Tables  1  and  2 
Cell  voltage  (V)  Grid 


60  x  30  x  12 

80  x  40  x  16 

90x48  x  18 

Current  density 
(A  cm-2) 

CPU  time  Relative  error 

(100%) 

Current  density 
(A  cm-2) 

CPU  time  Relative  error 

(100%) 

Current  density 
(A  cm-2) 

CPU  time 

0.47 

0.866 

17324 

3.12 

0.885 

43989 

0.99 

0.894 

105114 

0.57 

0.645 

14657 

3.84 

0.661 

37643 

1.43 

0.670 

89178 

0.67 

0.433 

6113 

3.81 

0.444 

31270 

1.44 

0.450 

79827 

0.77 

0.239 

4688 

3.51 

0.245 

24386 

1.33 

0.248 

56722 

0.87 

0.075 

3467 

2.67 

0.077 

20043 

1.00 

0.077 

47062 
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Table  4 


Comparison  of  present  study  with  results  from  literature  based  on  the  parameters  in  Ref.  [17] 


Cell  voltage  (V) 

Method 

Numerical  [17] 

Present  study 

Experiment  [17] 

Current  density  (A  cm-2) 

Relative  error  (100%) 

Current  density  (A  cm-2) 

Relative  error  (100%) 

Current  density  (A  cm-2) 

0.47 

0.998 

8.71 

0.889 

3.16 

0.918 

0.57 

0.856 

23.70 

0.659 

4.77 

0.692 

0.67 

0.462 

4.74 

0.428 

11.75 

0.485 

0.77 

0.179 

26.03 

0.219 

9.50 

0.242 

0.87 

0.035 

12.90 

0.032 

3.23 

0.031 

method.  The  influences  of  channel  aspect  ratio  and  GDL  thick¬ 
ness  are  described  in  detail  in  the  following  discussion. 

3.1.  Effects  on  cell  performance 

Flow  channel  geometry  is  important  for  optimizing  cell  per¬ 
formance  with  delivery  of  sufficient  reactant  to  the  reaction  sites. 
Various  configurations  with  different  channel  aspect  ratios  are 
investigated  to  explore  the  effects  of  channel  geometry  design. 
The  aspect  ratio  (AR)  is  defined  as  the  height  (in  Z-direction)  of 
the  channel  divided  by  the  width  (in  y-direction)  of  the  channel. 
Channel  height  and  width  are  modified  according  to  the  AR, 
such  that  the  channel  cross-sectional  area  is  fixed.  Hence,  the 
reactant  mass  flow  rate  at  the  inlet  and  the  channel  hydraulic 
diameter  are  the  same  in  all  case  studies.  To  maintain  a  fixed 
cell  width,  the  shoulder  width  is  also  changed.  Detailed  chan¬ 
nel  geometries  for  each  considered  AR  are  listed  in  Table  5. 
Fig.  2(a)  shows  the  relationships  between  output  current  den¬ 
sities  and  AR  at  moderate  and  low  cell  voltages.  Clearly,  cell 
performance  dependence  on  the  channel  aspect  ratio  varies  with 
operating  conditions.  The  cell  output  current  density  is  larger  at 
moderate  cell  voltage  and  performance  declines  rapidly  in  the 
high  reaction  rate  region,  when  the  AR  is  large  as  in  the  case 
of  1.5.  However,  despite  poor  performance  at  cell  potentials  of 
0.42-0.62  V  when  the  channel  is  flat,  such  as  at  AR  =  0.5,  current 
density  exceeds  that  obtained  with  other  channel  geometries  at 
0.14  V. 

Fig.  2(b)  displays  cell  current  density  curves  with  five  GDL 
thicknesses.  A  thinner  diffusion  layer  is  generally  believed  to 
reduce  the  reactant  vertical  diffusion  path  such  that  more  oxy¬ 
gen  arrives  at  the  catalyst  layer  to  gain  more  current.  In  contrast, 
cell  performance  varies  with  increased  thickness  of  the  diffu¬ 
sion  layer  because  the  reactant  is  transported  through  a  longer 
distance  and  the  flooding  problem  occurs.  However,  the  data  in 
Fig.  2(b)  reveals  that  optimal  thickness  increases  from  152  pan  at 


Table  5 

Channel  and  shoulder  geometries  for  various  ARs  (unit:  mm) 


AR 

Channel  height 

Half  channel  width 

Half  shoulder  width 

0.5 

0.57 

0.57 

0.23 

0.75 

0.69 

0.46 

0.34 

1 

0.080 

0.40 

0.40 

1.25 

0.089 

0.36 

0.45 

1.5 

0.98 

0.33 

0.47 

0.62  V  to  254  pm  at  0.22  V,  beyond  which  potential,  it  decreases. 
The  tendency  is  slightly  different  from  the  work  of  Jeng  et  al.  [9] 
with  a  monotonously  decreased  performance  when  GDL  thick¬ 
ness  increases  for  a  low-porosity  GDL.  By  using  a  high-porosity 
GDL,  however,  the  trends  coincide. 

Fuel  cell  output  is  controlled  by  three  degrading 
mechanisms — loss  associated  with  reactant  activation  energy  at 
the  electrodes;  loss  associated  with  transport  of  the  current  in 
the  electrically  conducting  media,  and  loss  associated  with  the 
transport  of  reactant  gases  in  the  cell.  These  are  called  activation 
overpotential,  ohmic  overpotential,  and  concentration  overpo¬ 
tential,  respectively.  The  AR  and  GDL  thickness  influences 
relative  magnitudes  of  these  three  irreversibilities  according  to 
the  operating  voltage,  as  is  obvious  from  the  data  in  Fig.  2. 
Local  variations  of  these  factors  may  also  play  an  important 


Fig.  2.  Current  density  distribution  at  various  cell  voltages  as  function  of  (a) 
AR;  (b)  GDL  thickness. 
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role  on  the  trend  transition  shown  in  the  figure.  A  fixed  cell 
overpotential  is  notably  designated  on  the  cathode  BP  surface 
in  the  present  study.  The  current  conservation  equation  solution 
provides  detailed  information  of  potential  distribution  through¬ 
out  the  domain.  Therefore,  ohmic  overpotential  and  activation 
overpotential  can  be  resolved  and  compared  with  the  oxygen 
concentration  for  different  operating  voltage.  Without  general¬ 
ity  loss,  representative  locations  in  the  cell  domain  are  assigned 
to  investigate  the  effects  of  these  two  parameters  and  various 
transport  properties  on  cell  performance  variation  in  the  follow¬ 
ing  sections. 

3.2.  Effects  on  membrane  conductivity  and  phase  potential 

A  series  of  demonstrations  involving  the  most  important 
transport  quantities,  such  as  phase  potentials,  reactant  concen¬ 
trations  at  selected  locations,  and  cell  voltages  are  performed  to 
elucidate  the  cause  of  the  forgoing  results.  All  of  the  data  are 
mirrored  to  yield  results  for  a  complete  channel  and  a  pair  of 
half  shoulders  adjacent  to  the  channel.  Notably,  the  negatively 
charged  electrons  move  from  low  potential  to  high  potential  and 
the  positively  charged  protons  move  oppositely.  The  main  con¬ 
cern  in  the  potential  variation  discussion  is  that  the  passages  with 
lowest  electron  potential  increase  from  the  cathode  boundary, 
and  the  ionic  potential  decrease  from  the  anode  catalyst  layer. 
Therefore,  ohmic  overpotential  decreases,  and  the  absolute  value 
of  the  electrochemical  reaction  driving  force,  the  activation  over¬ 
potential,  in  Eq.  (14)  increases.  Fig.  3  plots  and  compares  the 


(b)  Y  (m) 


Fig.  3.  Distributions  of  potential  and  conductivity  of  membrane  in  the  transverse 
direction  of  the  middle  X-Y  plane  for  three  values  of  AR  at  cell  voltages  of  (a) 
0.62  V;  (b)  0.14  V.  ((A)  AR  =  0.5,  (O)  AR  =  1  and  (□)  AR  =  1 .5).  The  short  bars 
are  the  interfaces  between  channel  and  shoulder. 


membrane  potential  and  conductivity  in  the  transverse  direc¬ 
tion  of  the  middle  X-Y  plane  of  the  membrane  for  AR  =  0.5,  1, 
and  1.5,  and  cell  voltages  of  0.62  and  0.14  V  at  Z- coordinate 
of  2.152  mm.  The  short  bars  on  the  plot  indicate  the  interfaces 
between  channel  and  shoulder  for  each  AR.  At  a  cell  voltage  of 
0.62  V,  the  membrane  conductivity  patterns  in  the  channel  region 
and  the  shoulder  differ,  according  to  the  AR  value.  At  the  mem¬ 
brane  location  of  the  channel  region,  conductivity  is  greater  at  a 
lower  AR.  Membrane  conductivity  is  a  function  of  water  activ¬ 
ity,  so  increasing  the  channel  transverse  dimension  at  AR  =  0.5 
facilitates  water  transport  at  the  anode  through  the  channel  to  the 
membrane,  causing  water  activity  in  the  central  region  to  exceed 
that  in  the  shoulder  region.  Trends  are  reversed  at  the  other  AR 
values  because  of  local  current  density  effect.  Membrane  phase 
potential  variations  are  strongly  related  to  the  conductivities;  a 
lower  conductivity  is  responsible  for  a  larger  membrane  ohmic 
loss  at  a  moderate  reaction  rate  at  which  the  local  change  in  cur¬ 
rent  density  is  expected  to  be  small.  However,  the  data  in  Fig.  3(b) 
indicate  that  at  a  high  reaction  rate,  higher  membrane  conduc¬ 
tivity  locations  exhibit  a  larger  membrane  ohmic  loss  because 
the  local  current  density  varies  markedly  at  this  cell  voltage. 
The  ohmic  law  and  the  fact  that  a  large  current  density  variation 
outweighs  the  trivial  local  conductivity  fluctuation  demonstrate 
that  ohmic  loss  in  the  membrane  phase  is  consistent  with  local 
current  density. 

Fig.  4  plots  the  membrane  phase  conductivity  and  potential 
at  GDF  thicknesses  of  152,  203  and  356  [xm  and  cell  voltages 
of  0.62  and  0.14  V,  to  examine  GDF  thickness  and  cell  volt- 


(b)  Y  (m) 


Fig.  4.  Distributions  of  potential  and  conductivity  of  membrane  in  the  transverse 
direction  of  the  middle  X-Y  plane  for  three  values  of  GDL  thickness  at  cell 
voltage  of  (a)  0.62  V;  (b)  0.14  V.  ((A),  (O)  and  (□)  are  denoted  as  152,  203  and 
356  pan). 
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age  effects  on  the  various  transport  properties.  The  data  were 
obtained  from  the  same  position  as  in  Fig.  3.  A  thick  GDL  sta¬ 
bilizes  transverse  variation  of  local  conductivity.  Reducing  this 
parameter  leads  to  a  highly  non-homogeneous  local  conductiv¬ 
ity  distribution  at  a  cell  voltage  of  0.14  V.  The  data  in  Fig.  4 
indicates  that  conductivity  is  related  to  membrane  potential  in  a 
manner  similar  to  that  in  Fig.  3.  At  a  low  reaction  rate,  these  two 
variables  vary  oppositely,  but  at  a  high  reaction  rate,  they  vary 
similarly. 

3.3.  Effects  on  solid  phase  potential  and  activation 
overpotential 

Solid  phase  potential  and  activation  overpotential  in  the  trans¬ 
verse  direction  for  various  AR  at  the  interface  between  the 
cathode  catalyst  layer  and  the  GDL  are  described  in  Fig.  5.  The 
data  in  Fig.  5(a)  reveal  that  the  shoulder  area  exhibits  a  small 
ohmic  overpotential  from  the  outer  surface  of  the  cathode  BR 
The  longer  electron  passage  in  the  catalyst  layer  beneath  the 
channel  region  corresponds  to  greater  potential  variation  on  the 
way  to  this  location.  Consequently,  the  activation  overpotential 
absolute  value  in  the  shoulder  region  exceeds  that  in  the  channel 
region,  indicating  that  the  electrochemical  reaction  driving  force 
is  stronger  there.  This  finding  of  non-uniform  activation  over¬ 
potential  is  consistent  with  that  of  Kulikovsky  et  al.  [8]  and  Sun 
et  al.  [11].  Comparing  the  potential  variation  effects  in  the  three 


selected  designs  of  the  channel  AR,  show  that  a  slender  chan¬ 
nel  design  exhibits  a  smaller  ohmic  overpotential,  and  a  larger 
absolute  activation  overpotential  at  all  positions  of  interest.  This 
result  is  explained  by  the  wide  electron  transport  passage  along 
the  shoulder  height  and  the  reduced  average  distance  between 
the  shoulders  and  the  channel  center,  leading  to  a  reduced  current 
resistance.  Variation  of  these  two  potentials  at  0. 14  V  appears  ini¬ 
tially  to  be  similar  to  that  in  Fig.  5(a).  The  channel  with  AR  =1.5 
has  a  stronger  driving  force  of  electrochemical  reaction  in  the 
shoulder  region  as  it  has  the  smallest  membrane  and  solid  phase 
ohmic  overpotentials  compared  to  the  other  two  channel  designs. 

Fig.  6  compares  the  transverse  variation  in  the  solid  phase 
potential  and  the  activation  overpotential  at  the  location  between 
the  cathode  GDL  and  the  catalyst  layer  at  three  GDL  thick¬ 
nesses  and  two  operating  potentials.  The  GDL  thickness  effects 
on  ohmic  overpotential  is  clearly  demonstrated  in  the  data  in 
Fig.  6(a).  Despite  the  fact  that  the  catalyst  layer  shoulder  area 
exhibits  minor  ohmic  overpotential  in  the  design  with  a  thickness 
of  152  pan,  the  potential  increases  abruptly  toward  the  channel 
central  region,  because  the  height  and  cross-sectional  area  of  the 
electron  transverse  transport  passage  are  small  in  this  thinnest 
design.  For  the  design  with  the  thicker  GDL,  the  larger  ohmic 
overpotential  in  the  catalyst  layer  shoulder  region  is  problematic, 
but  the  moderate  potential  variation  toward  the  central  chan¬ 
nel  is  advantageous.  Accordingly,  the  activation  overpotential 
variation  exhibits  the  same  tendency  as  that  of  the  solid  phase 
potential.  These  situations  are  similar  at  the  two  operating  volt- 


Fig.  5.  Transverse  distributions  of  solid  phase  potential  and  activation  over¬ 
potential  at  the  interface  between  the  cathode  catalyst  layer  and  the  GDL  for 
three  values  of  AR  at  cell  voltage  of  (a)  0.62  V;  (b)  0.14  V.  ((A)  AR  =  0.5,  (O) 
AR  =  1  and  (□)  AR  =  1 .5).  The  short  bars  are  the  interfaces  between  channel  and 
shoulder. 


Fig.  6.  Transverse  distributions  of  solid  phase  potential  and  activation  overpo¬ 
tential  at  the  interface  between  the  cathode  catalyst  layer  and  the  GDL  for  three 
GDL  thicknesses  at  cell  voltage  of  (a)  0.62  V;  (b)  0.14  V.  ((A),  (O)  and  (□)  are 
denoted  as  152,  203  and  356  |xm). 
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age,  but  the  absolute  values  and  variations  are  greater  in  the 
0.14  V  case. 

3.4.  Effects  on  water  saturation  level 

Channel  geometry  effects  on  the  transverse  saturation  level 
of  the  cathode  GDL  along  the  channel  direction  is  presented  in 
Fig.  7.  The  flat  channel  advantage  is  evidenced  by  liquid  water 
accumulation.  The  relatively  short  distance  through  which  water 
is  transported  at  AR  =  0.5  leads  to  the  low  saturation  level  in  the 
channel  direction  even  at  a  high  reaction  rate.  That  means  the 
fast  diffusion  associated  with  the  flat  channel  design  causes  the 
local  partial  pressure  of  water  vapor  to  be  low.  In  contrast,  sat¬ 
uration  level  at  the  rear  section  of  the  channel  with  AR  =1.5  is 
high,  as  the  water  vapor  cannot  easily  escape  from  GDL  under 
the  shoulder.  In  this  scenario,  effective  pore  space  in  the  porous 
medium  is  reduced  considerably,  and  more  mass  transport  over¬ 
potential  is  activated,  so  cell  performance  is  drastically  degraded 
at  low  cell  voltage. 

Fig.  8  plots  the  saturation  level  in  the  transverse  direction  of 
the  cathode  GDL  at  three  positions  v  =  0.05,  0.025  and  0.045  m, 
to  elucidate  GDL  thickness  effects  on  liquid  water  distribution. 
The  plot  demonstrates  that  the  thinner  GDL  design  has  the  lowest 
saturation  level  since  water  vapor  produced  in  the  cell  reaction 
easily  transports  to  the  channel  and  outside  the  cell.  On  the  con¬ 
trary,  the  longer  path  and  water  vapor  lower  diffusion  rate  in 
the  design  with  the  356  [xm- thick  GDL  results  in  greater  water 
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Fig.  7.  Effect  of  channel  aspect  ratio  on  transverse  water  saturation  of  the  cath¬ 
ode  GDL  along  channel  direction  at  cell  voltage  of  (a)  0.62  V;  (b)  0.14  V.  ((A) 
AR  =  0.5,  (O)  AR  =  1  and  (□)  AR  =1.5).  Right  half  portion  is  the  enlargement 
of  the  crowded  lines. 
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Fig.  8.  Effect  of  GDL  thickness  on  transverse  water  saturation  of  the  cathode 
GDL  along  channel  direction  at  cell  voltage  of  (a)  0.62  V;  (b)  0.14  V.  ((A), 
(O)  and  (□)  are  denoted  as  152,  203  and  356  |xm).  Right  half  portion  is  the 
enlargement  of  the  crowded  lines. 


vapor  concentration  with  a  high  probability  of  over- saturation 
and  liquid  water  formation,  contributing  to  concentration  over¬ 
potential  and  performance  reduction,  especially  at  v  =  0.45  m  at 
a  cell  potential  of  0. 14  V. 

Note  that  temperature  variation  is  not  considered  in  this  study, 
model  result  may  over  predict  saturation  level  at  high  reaction 
rate  because  the  saturation  pressure  would  be  greater  accom¬ 
panied  by  larger  local  temperature.  However,  the  membrane 
water  content  adjacent  to  the  high  temperature  region  is  ele¬ 
vated  due  to  isothermal  assumption.  This  could  enhance  ionic 
current  transport  and  compensate  the  loss  from  over-predicted 
saturation  level.  Nevertheless,  this  highlights  the  necessity  of 
conducting  a  non-isothermal  analysis  to  interpret  these  mutual- 
coupled  effects. 


3.5.  Effects  on  oxygen  concentration  and  local  current 
density 

Fig.  9  depicts  oxygen  mass  fraction  distribution  and  local  cur¬ 
rent  density  under  the  same  condition  and  at  the  same  location  as 
in  Fig.  5.  Intuitively,  the  flat  channel  supplies  more  oxygen  to  the 
catalyst  layer  and  a  high  current  density  is  expected.  However, 
calculations  reveal  entirely  different  trends  between  channel  and 
shoulder  region  catalyst  layers.  The  design  with  AR  =  0.5  gen¬ 
erates  the  lowest  local  current  density  in  the  channel  region, 
despite  its  having  the  largest  oxygen  mass  fraction,  and  the 
expected  relationship  between  reactant  concentration  and  the 
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Fig.  9.  Transverse  distributions  of  oxygen  mass  fraction  and  local  current  density 
at  the  interface  between  the  cathode  catalyst  layer  and  the  GDL  for  three  values 
of  AR  at  cell  voltage  of  (a)  0.62  V;  (b)  0.14  V.  ((A)  AR  =  0.5,  (O)  AR=  1  and 
(□)  AR=  1.5). 


Fig.  10.  Transverse  distributions  of  oxygen  mass  fraction  and  current  density 
at  the  interface  between  the  cathode  catalyst  layer  and  the  GDL  for  three  GDL 
thicknesses  at  cell  voltage  of  (a)  0.62  V;  (b)  0. 14  V.  ((A),  (O)  and  (□)  are  denoted 
as  152,  203  and  356  |xm). 


reaction  rate  is  not  observed  except  at  the  shoulder  center.  On 
average,  the  slender  design  with  AR  =  1.5  generates  more  cur¬ 
rent  at  moderate  cell  voltage.  This  phenomenon  is  explained 
by  the  fact  that  in  this  scenario,  the  reaction  is  relatively  slow 
and  the  high  reactant  concentration  is  not  as  important  as  low 
ohmic  overpotential  and  high  activation  overpotential  provided 
by  slender  channel  geometry  such  as  AR  =  1.5.  This  design  has 
a  wider  rib  zone  than  those  of  other  designs,  providing  a  small 
increase  and  variation  in  the  solid  phase  potential  as  well  as  large 
absolute  activation  overpotential,  causing  high  current  density 
and  favorable  cell  performance.  At  a  high  reaction  rate  as  shown 
in  Fig.  9(b),  the  dominant  mechanism  of  local  current  density 
transits  from  electric  potential  at  the  channel  region,  to  oxygen 
concentration  at  the  shoulder  region  at  AR=  1.5.  This  transition 
is  also  found  when  AR  =  1 .  The  design  with  AR  =  0.5  exhibits  the 
same  local  current  density  variation  trend  as  in  the  case  of  0.62  V. 
At  AR=  1.5,  a  high  reactant  concentration  need  outweighs  a 
high  activation  overpotential  need,  so  the  expected  relationship 
between  the  concentration  and  the  local  reaction  rate  appears 
earlier.  The  large  shoulder  area  width  hinders  reactant  transport, 
resulting  in  a  relatively  low  level  of  local  oxygen  mass  fraction 
and  a  sharp  decline  in  the  local  current  density  under  the  shoul¬ 
der  area.  In  contrast,  the  sufficient  oxygen  provided  by  the  flat 
channel  such  as  AR  =  0.5  causes  most  of  the  region  to  exhibit 
a  potential  controlled  state  and  on  average,  produces  a  greater 
current  than  that  generated  by  other  designs. 

GDL  thickness  effects  on  transverse  distribution  of  oxygen 
concentration  and  current  density  is  plotted  in  Fig.  10.  Exhibiting 


a  trend  opposite  to  that  of  the  potential,  a  thinner  GDL  provides 
more  oxygen  to  the  channel  region  catalyst  layer.  The  thickest 
design  is  associated  with  a  greater  concentration  in  the  shoul¬ 
der  region.  These  variations  are  related  to  the  vertical  depth 
and  transverse  cross-section  of  reactant  delivery.  In  the  case 
in  which  GDL  thickness  is  152  pan,  the  vertical  path  is  short 
and  oxygen  concentration  is  high  at  the  catalyst  layer  beneath 
the  channel.  Nevertheless,  the  transverse  transport  cross-section 
is  reduced  and  oxygen  concentration  falls  substantially  at  the 
shoulder  region.  With  reference  to  local  current  density  dis¬ 
tribution,  mechanisms  of  activation  overpotential  mechanism 
and  oxygen  concentration  have  different  effects  in  different 
regions.  The  activation  overpotential  dominates  the  reaction  in 
the  catalyst  layer  under  the  channel  region  because  oxygen  con¬ 
centration  fulfills  the  requirement  for  electrochemical  reaction. 
Accordingly,  the  current  density  trend  is  consistent  with  that  of 
the  activation  overpotential  throughout  the  entire  region  in  the 
356  jum  case  and  in  most  of  the  region  in  the  203  and  152  pm 
cases  at  0.62  V.  A  peak  point  in  the  concave  pattern  of  the  local 
current  density  appears  only  near  the  central  shoulder  region  in 
the  152  pm  cases,  suggesting  that  from  this  point  to  the  central 
shoulder,  oxygen  deficiency  forces  the  current  density  to  drop 
according  to  the  oxygen  mass  fraction,  and  local  performance  is 
dominated  by  oxygen  concentration.  Nevertheless,  most  current 
density  gain  in  the  shoulder  region  catalyst  layer  arises  from 
higher  activation  overpotential  for  the  design  with  the  152  pm- 
thick  GDL,  compensating  the  loss  at  the  central  region  and 
generating  more  current  than  other  designs  at  0.62  V. 
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Similar  result  can  be  found  in  the  work  of  Sun  et  al.  [  1 0]  which 
investigated  the  GDL  thickness  effects  at  total  call  overpotentials 
of  —0.4,  —0.5  and  —0.6  V.  The  thinner  GDL  thickness  design 
exhibits  the  slowest  reaction  rate  at  the  catalyst  site  under  the 
channel  region.  This  phenomenon  was  interpreted  by  a  longer 
electron  traveling  length  to  the  channel  region  than  to  the  shoul¬ 
der  region,  such  that  fewer  electrons  would  participate  in  the 
electrode  reaction  in  the  region  under  the  channel. 

At  a  high  reaction  rate  with  a  cell  voltage  of  0. 14  V  as  shown 
in  Fig.  10(b),  oxygen  concentration  is  lower  than  in  Fig.  10(a) 
because  the  reaction  rate  is  higher,  but  transport  depth  and 
transverse  cross-section  effects  on  oxygen  concentration  are  the 
same.  However,  local  current  density  distributions  exhibit  quite 
different  patterns.  The  plot  depicts  clear  peak  points  in  the  con¬ 
cave  pattern  for  all  three  GDL  thicknesses.  As  stated  previously, 
the  356  jam  design  generates  more  current  than  the  other  two 
designs  at  the  central  channel  region  catalyst  layer  because  the 
activation  overpotential  dominates  the  reaction.  This  situation 
changes  rapidly  toward  the  shoulder  region  since  oxygen  con¬ 
centration  is  extremely  low.  The  figure  also  shows  that  the  peak 
current  density  points  move  near  the  shoulder  in  the  remain¬ 
ing  two  cases.  A  thinner  GDL  corresponds  to  a  latter  peak 
point  from  the  channel  central  and  a  larger  peak  current  den¬ 
sity.  Beyond  these  points,  current  density  variation  is  nearly  the 
same  as  oxygen  concentration.  At  the  peak  points,  the  design 
with  the  152[xm-thick  GDL  clearly  has  the  largest  local  cur¬ 
rent  density;  unfortunately,  oxygen  concentration  falls  sharply 
toward  the  shoulder  center  and  current  density  is  the  lowest  of 
the  three  designs.  Consequently,  the  152  ptm-thick  case  does  not 
exhibit  the  best  cell  performance.  On  the  other  hand,  despite  the 
fact  that  the  203  jxm-thick  design  has  a  moderate  local  current 
density  output  at  the  channel  center,  peak  point,  and  shoulder 
center,  its  local  current  density  curve  envelopes  the  largest  area 
with  the  X—Y  plane  and  becomes  the  optimal  GDL  thickness  at 
low  cell  voltage.  This  finding  is  quite  different  from  the  result 
of  Sun  et  al.  [10]  that  no  significant  variation  of  average  current 
density  was  found  when  GDL  thickness  reduced  from  0.35  to 
0.15  mm.  The  discrepancy  may  be  attributed  to  different  oper¬ 
ating  conditions  or  cell  geometries  between  these  two  models. 

4.  Conclusions 

A  multi-dimensional,  multi-component,  computational  fluid 
dynamic  model  was  employed  to  elucidate  transport  phenom¬ 
ena  and  electrochemical  reaction  in  PEMFCs.  Two  important 
transport  component  design  parameters — channel  aspect  ratio 
and  GDL  thickness,  are  investigated  in  detail.  The  simulation 
and  discussion  support  the  following  conclusions: 

1 .  The  designs  of  channel  aspect  ratio  and  GDL  thickness  affect 
the  physical  property  distributions.  At  a  high  reaction  rate, 
these  two  parameters  have  a  strong  influence  on  cell  perfor¬ 
mance. 

2.  Three  mechanisms  of  cell  irreversibility  are  resolved  locally 
from  the  model  and  are  considered  to  determine  variation  in 
the  macroscopic  cell  current  density  and  performance. 


3.  Transverse  current  density  distribution  is  governed  by  both 
electron  conduction  and  activation  overpotential  in  channel 
region  or  by  reactant  concentration  in  shoulder  region  of  cat¬ 
alyst  layer.  The  relative  strengths  of  these  two  factors  depend 
on  the  operating  voltage  and  the  transport  geometry. 

4.  At  a  moderate  reaction  rate,  the  transverse  direction  current 
density  in  most  regions  is  a  convex  function  of  position  and 
is  influenced  by  solid  phase  potential  and  activation  overpo¬ 
tential,  making  the  peak  point  of  the  concave  pattern  current 
density  either  close  to  the  shoulder  region  or  non-existent. 
Thus,  a  geometry  design  that  facilitates  electron  transport 
such  as  a  large  channel  aspect  ratio  or  a  thin  GDL  thickness 
causes  a  relatively  larger  current  in  the  catalyst  layer. 

5 .  At  a  low  cell  voltage,  the  largest  reaction  rate  location  is  close 
to  the  channel  center,  therefore,  the  electrochemical  reaction 
in  the  majority  regions  is  dominated  by  the  reactant  transport 
and  a  smaller  channel  aspect  ratio  design  is  preferred. 

6.  Large  GDL  thickness  positively  increases  oxygen  transport 
under  the  shoulder  region;  therefore,  output  current  density 
elevates  according  to  cell  voltage  decrease.  However,  at  the 
lowest  considered  cell  voltage  of  0.14  V,  oxygen  deficiency 
caused  by  long  traveling  length  and  liquid  water  clogging 
effect  reverses  this  relationship. 
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